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Localized Wannier functions provide an efficient and intuitive means by which to compute dielec- 
tric properties from first principles. They are most commonly constructed in a post-processing step, 
following total-energy minimization. Nonorthogonal generalized Wannier functions (NGWFs) 1 may 
also be optimized in situ, in the process of solving for the ground-state density. We explore the 
relationship between NGWFs and orthonormal, maximally localized Wannier functions (MLWFs) 2 , 
demonstrating that NGWFs may be used to compute electric dipole polarizabilities efficiently, with 
no necessity for post-processing optimization, and with an accuracy comparable to MLWFs. 
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In this Brief Report, we explore the equivalence 
between nonorthogonal generalized Wannier functions 
(NGWFs) 1 , generated using linear-scaling Kohn-Sham 
density functional theory (DFT) 3 , and their orthonormal 
counterparts, particularly maximally localized Wannier 
functions (MLWFs) 2 , both recently reviewed in Ref. 4. 
We demonstrate the comparable, high accuracy of the 
two formalisms for dielectric response, laying the foun- 
dation for large-scale calculation of optical properties. 

We begin with the single-particle density-matrix de- 
fined, for a given set of Bloch orbitals \ip n k), where n 
indexes occupied bands, k is the crystal wave- vector and 
we suppress the spin index for notational clarity, by 
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Here, 1BZ is the first Brillouin zone corresponding to 
the periodic unit cell of volume V cc \\. A reformulation 
of such Bloch states suitable for the study of spatially 
localized properties was proposed by Wannier 5 , whose 
eponymously named functions are defined, for a unit cell 
at the lattice vector R, by 
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The orthonormality of Bloch orbitals is preserved, 

(w„R,|m m R/) = 5 nm 5 RR >, 
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and we may choose the gauge freely, so that any prior 
unitary transformation among the orbitals, \ip n k) = 
SmlV , mk)t/mnk, may also give rise to a valid set of gen- 
eralized Wannier functions, via Eq. 2. Unoccupied states 
may be included in the wannierization, while maintaining 
the same occupied density, by appropriately transform- 
ing the occupancy of the orbitals to give / nk , to give 



/«mk = J2 P Ul pk fpkUpmk- The density-matrix may be 
readily expressed in terms of Wannier functions, in the 
separable form proposed in Ref. 6, and given by 

P = y^JwnR)knmR'-R(w m R'\, where (4a) 
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is commonly known as the density kernel. 

The extension of this formalism to nonorthogonal gen- 
eralized Wannier functions (NGWFs) is both of practi- 
cal interest and utility. Orthonormality and spatial lo- 
calization are generally competing requirements 7 , hence 
nonorthogonal orbitals may form a more efficient basis 
in which to expand short-ranged operators and, as a 
result, they are used extensively in linear-scaling DFT 
approaches. We may express these NGWFs, \<j) a n), 
simply in terms of the generalization of the transfor- 
mation matrices t/ k to possible non-unitarity matrices 
M k , that is M na u = ("0nk|^ak), whereafter / Q/3k = 
J2 n M^ nk f n kM n pk. We use Latin and Greek letters to 
index orthonormal and nonorthogonal sets, respectively, 
and implicitly sum over repeated index pairs. 

In the nonorthogonal case, the density-matrix may be 
expanded in separable form via the tensor contraction 
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and the price to be paid for nonorthogonality is a nontriv- 
ial metric tensor given by S a /3 — (</>aR|<A3R')6RR', which 
defines the inter-relationship between covariant vectors, 
|<^aR) = \4>n)Si3a, and contravariant vectors (NGWF du- 
als) \4>^) = \4>/3R)SP a . The contravariant metric S a ^ in 



Eq. 5b is defined such that 5 , Q7 S' 7 ' 3 = 8J 3 , and is also in- 
dependent of the lattice vector. Orthonormality is thus 
replaced by the general tensor expression 

<0aR|</> 7 R<>S 7/3 = S a7 (^|<,) = 5jS RR ,. (6) 

Numerous optimization procedures have been devel- 
oped for ab initio Wannier functions. A widespread ap- 
proach involves their construction in a post-processing 
step, computing the t/„ m k or M na \^ matrices, and then 
the generalized occupancies /k, following the computa- 
tion of the delocalized orbitals. However, it has been rec- 
ognized, and utilized in the context of large-scale calcu- 
lations for some time 8 , that localized Wannier functions 
may also be optimized directly in situ, that is during 
the process of solving for the electronic structure. In the 
latter, the basis expansion of the functions, {0 a R.( r )} 
and the corresponding density kernel K R must be op- 
timized together, reconstructing the delocalized orbitals 
afterwards only if necessary. 

A variety of plausible criteria may also be employed for 
Wannier function optimization in either case, such as en- 
ergy downfolding 9 or maximal Coulomb repulsion 10 , or, 
as used in this work, total-energy minimization or spatial 
localization. Depending on their definition, these criteria 
may or may not uniquely define the Wannier functions, 
in that they may admit some residual gauge freedom. A 
particularly efficacious measure for localization is the sec- 
ond central moment which, for a set of Wannier functions 
{\w n -R.}} takes the form of the spread functional, 
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where the generalization to the nonorthogonal case does 
not yield straightforward physical interpretation. ML- 
WFs 2 are those orthonormal Wannier functions gener- 
ated by unitary transformations f/ ram k that minimize $7, 
for a fixed set of orbitals. MLWFs are usually computed 
in a post-processing procedure, using an implementation 
such as WANNIER90 11 , and have been widely adopted 
as an accurate minimal basis with which to compute nu- 
merous ground-state and excited-state properties, as well 
as to augment DFT with many-body interactions 4 . ML- 
WFs have been used to great effect, moreover, in the 
context of molecular dynamics, particularly interesting 
examples including the calculation of the dielectric per- 
mittivity and dipolar correlation of liquid water 12 , as well 
as its dynamical charge and dipole tensors 13 . 

NGWFs, unlike their orthonormal counterparts, are 
more commonly optimized in situ, as a by-product of 
total-energy minimization with respect to the density- 
matrix, for example, in the ONETEP linear-scaling 
DFT code 1 ' 14 . In the latter, NGWFs are expanded in 
a fixed underlying basis of periodic cardinal sine func- 
tions (also known as psinc 15 or band-width limited S- 
functions), whose spatial finesse is determined by a sin- 
gle variational parameter, the kinetic energy cutoff of the 
equivalent plane-wave basis. The NGWFs are then those 




FIG. 1. (Color online) Wannier functions, NGWF (left) and 
MLWF (right), of predominantly oxygen p z 2 (highest occu- 
pied, l&i) character in H2O, at zero applied field, with iso- 
surfaces at one sixth of their respective maxima. Both types 
retain some residual arbitrariness following optimization. 



functions, when traced with their corresponding opti- 
mized density kernel, which reproduce the ground-state 
density-matrix, whence the ground-state energy 
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In practice, in order to extremize the total-energy 
with respect to idempotent density matrices, two nested 
conjugate-gradients variational minimization procedures 
are performed. In the inner loop, the energy is minimized 
with respect to the elements of the density kernel, for a 
fixed NGWF expansion, and in the outer, the density 
kernel is kept fixed while the total energy is minimized 
with respect to the NGWF psinc-expansion. A number of 
similar methods have been proposed in which equations 
of motion generate optimized nonorthogonal functions 8 . 

An intuitive interpretation of Wannier functions is fur- 
nished via the modern theory of polarization 16 , in that 
changes in their centers {r) nm = {w n o\r\w m o) exactly re- 
produce, and thus may be used to efficiently calculate, 
changes in the polarization of insulating systems. The 
change in electronic polarization <SP, subject to a gap- 
preserving perturbation, may be expressed as 
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where the k-independence of the occupancies (also spin- 
degenerate) implies that it is sufficient to consider only 
the R = term. Here, respectively, we have provided the 
orthonormal case for N occupied bands, and the more 
general, nonorthogonal case. 

It has been shown that close-to-orthonormal Wannier 
functions generated by means of direct minimization, of 
an appropriately constructed functional, may be used to 
efficiently compute dielectric properties 17 . It is of impor- 
tance, particularly for linear-scaling methods, to general- 
ize this result and verify that in situ optimized NGWFs 



a 


NGWF 


MLWF 


Gaussian 


Experiment 


H 2 


10.58 


10.47 


10.76", 7.4° 


9.64°, 9.79-^ 


NH 3 


15.28 


15.24 


15.63°, 12.1° 


14.56°, 18. 9 e 


CH 4 


17.70 


17.49 


17.74°, 14.8° 


17.27°, 17.5 eJ 


C2H4 


28.50 


28.39 


28.77°, 25.6° 


27.70 d , 28.69 / 


CO 


13.64 


13.53 


13.73°, 12.1° 


13.09°, 12.8 e , 13.16 / 


CO2 


18.00 


17.85 


18.06°, 14.8° 


17.51°, 19.6 e - / , 17.48 / 


N 2 


11.99 


11.87 


12.31°, 10.8° 


11.74° J , 11. 5 e 


CioHg 


122.8 


123.0 


121. 76 c 


117. 4? > v , 118. 9 Z 


K 


NGWF 


MLWF 


Gaussian 


Experiment 


H 2 


0.16 


0.14 


0.14° 


0.67° 


NH3 


2.45 


2.64 


2.70° 


1.94° 


C2H4 


12.18 


12.03 


11.94° 


11.4° 


CO 


3.53 


3.54 


3.55° 


3.57° 


co 2 


13.88 


13.96 


13.70° 


13.83 d , 13.70 / 


N2 


4.50 


4.55 


4.83° 


4.59°, 4.45 / 


CioHg 


96.1 


94.8 


94.13 c 


86.9*, 79.0 H , 63.6 Z 



( a ' Static polarizability in a d-aug-cpVTZ basis 19 . 

(i>) Static PBE polarizability in 6-311++G(d,p) basis at 

B3LYP/6-311G** optimised geometries 20 . 
' c ' Static polarizability in Sadlej pVTZ basis 21 . 
( d > Compiled in Ref. 19, based on analysis of anisotropy data 22 . 
( e ' CRC Handbook 23 . W Extrapolated Raylcigh scattering 24 . 
W Anisotropic refraction 25 . ( y > Laser Stark spectroscopy 25,26 . 
("' Optical measurements at 632.8 nm in solution 27 . 

TABLE I. Isotropic (a) and anisotropic (k) polarizabilities 
(e 2 a 2 , Ha -1 ), from DFT using nonorthogonal (NGWF) and 
orthonormal (MLWF) Wannier functions. Previous Gaussian- 
basis calculations 28 , and experimental values are included. 



can reproduce electronic response properties with the 
same reliability as that of the well documented MLWFs, 
as NGWFs are increasingly being used in large-scale 
methods for spectral partitioning and dielectric proper- 
ties, particularly in molecular systems 18 . The simplest 
such response property is perhaps the high-frequency 
(termed "clamped- ion" or "static" ) linear dipole polariz- 
ability tensor 
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where £ is an applied electric field within the dipole ap- 
proximation. This polarizability is somewhat different 
from that which is most frequently probed experimen- 
tally, namely the static or visual frequency regimes, and 
neglects the response of the ionic positions. 

Two different Kohn-Sham DFT packages were used in 
order to compute polarizabilities within the NGWF and 
MLWF formalisms, respectively the ONETEP linear- 
scaling code 1 ' , and a combination of a plane- wave pseu- 
dopotential package 29 and the WANNIER90 11 code. An 
example of each type of function is depicted in Fig.l. A 
set of well-isolated, closed-shell molecules were selected, 



so that a sawtooth-potential representation of the electric 
field could be used, with the potential boundary max- 
imally distant from the molecules, up to a maximum 



field value of ±8.0 x 10~ 5 Ha e" 



in intervals of 



2.0 x 10~ 5 Ha e _1 a^ , for all systems. The response re- 
mained well within the linear regime at these field values, 
which lay well below the threshold for Zener breakdown. 
The rates of change in polarization was calculated us- 
ing linear-regression of finite-difference data. Identical 
norm-conserving pscudopotcntials 30 were used in both 
cases, having been generated in the required formats, the 
Perdew-Burke-Ernzcrhof (PBE) exchange-correlation ap- 
proximation 28 , and the same run-time parameters and 
analysis were used for both codes so far as possible. 
Zero-field ground-state geometries were optimized using 
ONETEP, in cubic simulation cells of side length 40 ao 
(50 ao in the case of naphthalene CioHg). The density 
and potential were fully reset to those of atomic superpo- 
sitions upon each incrementation of the field. An equiva- 
lent plane- wave cutoff of 1000 cV, T-point Brillouin zone 
sampling, no density-kernel truncation and NGWFs with 
a 10 ao radius cutoff were used. 

In the case of nonaxially symmetric molecules, random 
initial guesses for the MLWFs were regenerated at each 
incrementation of the electric field. For axially symmet- 
ric molecules such as CO, CO2 and N2, however, the 
maximal localization condition does not uniquely define 
the MLWF centers under rotations about the axis, as dis- 
cussed in Ref. 31. While the sum of centers, and hence 
the transverse response, should be well-defined, in prac- 
tice this unbroken symmetry results in excessively noisy 
linear-response data. It was found, however, that re- 
initializing the MLWFs to s-orbitals at each field value, 
with centers coinciding with a set of zero-field MLWFs, 
proved sufficiently robust to obtain excellent linear fit- 
ting. No such measures were necessary in the case of 
ONETEP NGWFs, due to an effective symmetry break- 
ing introduced by the underlying real-space psinc grid. 

The isotropic and anisotropic parts of the polarizability 
tensor a are defined, respectively, as 
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our computed values of which using NGWFs and MLWFs 
are shown in Table I. The quadratic mean fractional dis- 
crepancy between the isotropic parts was 8. Ox 10~ 3 , while 
the discrepancy was greater for the anisotropic parts, at 
4.1 x 10~ 2 . As judged by the arithmetic mean frac- 
tional discrepancies (given henceforth in parentheses), 
the NGWFs tended to provide slightly larger isotropic 
parts (by 6.6 x 10~ 3 ), and also anisotropic parts (by 
1.3 x 10~ 3 ), than the MLWFs. Perhaps serendipitously, 
the NGWF values lay closer than the MLWF results, 
for the isotropic parts, in all cases, to the previous 
DFT(PBE) calculations of Ref. 19, computed using a 
sophisticated time-dependent coupled-perturbed method 
with a triple-^ Gaussian basis set; the quadratic (arith- 
metic) mean discrepancies with respect to these previ- 



cms results were 1.5 x 1(T 2 (1.3 x 1CT 2 ) and 2.2 x 1CT 2 
(2.0 x 1CT 2 ), respectively. The trend was reversed for 
anisotropies. 

Polarizabilities calculated using the related Wannicr 
function varieties agree rather well in spite of the the 
significant technical dissimilarities between the ab initio 
packages generating them, and there arc a number of pos- 
sible origins for the small discrepancies observed. First 
considering the NGWF and MLWF values, both based on 
the plane-wave formalism and using the same ionic ge- 
ometry, the NGWF method is the more approximate in 
that it spatially truncates the Wannier functions and the 
kinetic-energy operator. Moreover, these methods differ 
in their handling of pseudopotcntials, and, substantially, 
in their energy-minimization algorithms. With respect 
to the previous Gaussian-basis results of Rcf. 19, the 
possible origins for discrepancy are manifold, most no- 
tably, the ionic geometries employed differ and the latter 
method treats the core electrons explicitly. 

The probable errors (arising from the linear fit to the 
data) in the isotropic and anisotropic polarizabilities, de- 
noted Aa and An, respectively, and given by 
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were computed using the unbiased variance estimate 
(Actij) on each polarizability component a^, and are 
shown in Table II. The noise in the data for NGWFs is 
somewhat more system dependent, as the NGWF trun- 
cation depends on the ionic geometry, and higher than 
in the MLWF case for most of the molecules studied. 
The quadratic (arithmetic) mean of the ratio of the 
estimated error in the isotropic polarizability a to its 
value was estimated at 1 x 1CT 3 (4 x 1CT 4 ) for NGWFs 



and 2 x 1CT 4 (2 x 1CT 4 ) for MLWFs. Correspondingly, 
for the anisotropic part Aa, we estimated these ratios 
to be, respectively, 1 x 10~ 2 (7 x 10~ 3 ) and 3 x 1CP 3 
(2 x 10~ 3 ). Nonetheless, the probable errors in the linear 
fits to the polarizability data were extremely small using 
both methods, and inconsequential with respect to the 
expected errors in the approximate functional. 

In conclusion, we have shown that nonorthogonal Wan- 
nier functions optimized in situ may be used to com- 
pute molecular polarizabilities with an accuracy compa- 
rable to MLWFs post-processed from plane- wave DFT. 
This result is promising for the computation of numerous 
dielectric properties, and the full application of linear- 
scaling Wannier function analysis to large systems. A 
promising avenue for future work is the generalization of 
a method for the dielectric response in extended systems, 
such as that described in Ref. 32 and applied to solids in 
Ref. 17, to the linear-scaling NGWF formalism. 
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TABLE II. Probable fractional errors in molecular electric 
polarizabilities computed using NGWFs and MLWFs. 
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